Time Delay Effect on the Love Dynamical Model 
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We investigate the effect of time delay on the dynamical model of love. The local stability analysis 
proves that the time delay on the return function can cause a Hopf bifurcation and a cyclic love 
dynamics. The condition for the occurrence of the Hopf bifurcation is also clarified. Through a 
numerical bifurcation analysis, we confirm the theoretical predictions on the Hopf bifurcation and 
obtain a universal bifurcation structure consisting of a supercritical Hopf bifurcation and a cascade 
of period-doubling bifurcations, i.e., a period doubling route to chaos. 

PACS numbers: 02.30.Ks, 05.45.-a, 89.65.-s 

Keywords: Love dynamics, Delay differential equations, Bifurcation analysis 



I. INTRODUCTION 

In a one page pioneering paper [l| and a book Q , Stro- 
gatz suggested a simple pedagogical model describing a 
love affair. His goal was to teach harmonic oscillation 
phenomena using "a topic that is already on the minds 
of many college students: the time-evolution of a love af- 
fair between two people" . Later, Rinaldi, Gragnani, and 
Feichtinger proposed more realistic mathematical models 
for love dynamics : They showed that dynamic phe- 
nomena in the field of social science can also be analyzed 
by using a modeling approach via ordinary differential 
equations. Their models explained that two individuals, 
who are completely indifferent to each other from the 
start, approach a plateau of love affair. They also showed 
that the coexistence of insecurity and synergism results 
in a cyclic dynamics of romantic feelings. Secure individ- 
uals react positively to their partner's love and are not 
afraid about their partner becoming emotionally close 
to them while non-secure ones react negatively to high 
involvement Synergic individuals are those who in- 
crease their reactions to their partner's appeal when they 
are in love Q. Thereafter, following a suggestion of Stro- 
gatz, Sprott investigated the dynamics of a love triangle, 
which produce a chaotic behavior 0, and suggested dy- 
namical models of happiness ||. Moreover, Wauer et al. 
studied the love dynamics for time-varying fluctuations 
Q. Very recently, Rinaldi et al. constructed full catalog 
of possible love stories among two individuals [lfj, and 
Barley and Cherif studied the stochastic love dynamical 
model 

On the other hand, a dynamical model described by 
using delay differential equations (DDEs) has attracted 
much attention in various fields of science j_e.g-, biology 
[lH, [TJJ , chemistry [lij , neural syste ms fl5l fill ] , excitable 
systems fl7l . transport control [H, Il9l| and cryptogra- 
phy [2(| l2l| . DDEs also support a realistic mathemati- 
cal modeling of economic dynamics [22 24]. In the field 
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of social science, Liao and Ran recently showed that the 
time delay on love dynamics can cause a Hopf bifurcation 
25]. However, they did not consider the linear, secure, 
and non-secure returns. Also, they did not deal with the 
synergic instinct. For that reason, they could not clarify 
the condition for the occurrence of Hopf bifurcation. 

In this paper, we investigate the effect of time delay 
on the nonlinear dynamical model describing a love affair 
between two individuals. By analyzing the characteristic 
equation of linearization of the model, we theoretically 
prove that if no one exhibits a non-secure return in both 
cases of synergic and non-synergic couples, then the exis- 
tence of time delay cannot disturb a plateau of love affair, 
i.e., steady state. However, if at least one of them has a 
non-secure return, then the time delay on return function 
can cause a Hopf bifurcation and a cyclic love dynamics. 
On the other hand, through numerical bifurcation analy- 
sis, we confirm the theoretical results on Hopf bifurcation 
and investigate additional bifurcation phenomena. As a 
result, we obtain a universal bifurcation structure con- 
sisting of a supercritical Hopf and a cascade of period- 
doubling bifurcations, resulting in chaotic motion for our 
models. 

This paper is organized as follows: Section [TTI presents 
the love dynamical model. In Section IIII1 we prove the 
occurrence of Hopf bifurcation by using a local stability 
analysis. In Section IIV1 we show the results of our nu- 
merical bifurcation analysis. The conclusion is given in 
Section [Vj 



II. LOVE DYNAMICAL MODEL 

Now, let us investigate the love dynamical model based 
on a series of models proposed by Rinaldi, Gragnani, and 
Feichtinger The model has a variable Xi (i = 1,2), 

which is a measure of the love of an individual i for his 
or her partner j {J = 2,1). Positive values of Xi repre- 
sent love while negative values are associated with hate. 
Complete indifference is identified by Xi = 0. 

It is important to mention the time scale of the love 
dynamical model. Fast fluctuations of the feelings in- 
fluenced by daily or weekly activities cannot be captured 
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by the model. Also, the learning and the adaptation pro- 
cesses over a long range of time are not considered. Thus, 
the following model can only be used on an intermediate 
time scale (months/years), for example, in predicting if a 
love story will be characterized by stationary or stormy 
feelings [4]. 

The dynamics of love is comprised of three basic pro- 
cesses: oblivion O, , return i?, and instinct L. 



±,{1) = Oi(xi(t)) + R i (x j {t)) +Ii(xi(t)) 



(1) 



In the following, Xi(i) and Xj(t) are replaced by Xi and 
Xj for compact notation. Oblivion is described by 



(2) 



where > is a forgetting coefficient. Therefore, Xi 
decays exponentially when an individual i loses partner 
j(R i = I i = 0). 

The return R i is related to the reaction of individual i 
to the partner's love Xj. Three different types of return 
functions have been considered, namely, the linear return 
Ri , the secure return Ri 3 , and the non-secure return R4 1 . 
The simplest one is the linear return Q given by 



Ri (p^y) — PiXj 



(3) 



where > is a reactiveness to the love. It is unbounded 
and describes that an individual i "loves to be loved" and 
"hates to be hated" . The secure return Q is specified by 



j3iXj/{\ + Xj) for Xj > 0, 
/3iXj/(l — Xj) for Xj < 0. 



(4) 



It is an increasing and bounded function, as shown in Fig. 
Ufa). Secure individuals react positively to their part- 
ner's love and are not afraid about their partner becom- 
ing emotionally close to them. However, non-secure in- 
dividuals react negatively to high pressures and involve- 
ment, as shown in Fig. Hlb). The non-secure return Q 
is described by 



/^•(l-z/v^i + ^xi + x/)} 



fcXj/il-Xj) 



for x 3 > 0, (5) 
for Xj < 0. 



The instinct is related to the reaction of individual i 
to the partner's appeal Aj. In the following, we only con- 
sider positive appeal Aj. Two different types of instinct 
functions have been suggested, namely the synergic in- 
stinct Ii* and non-synergic instinct I™ . The non-synergic 
instinct [3, H| is given by 



(0) 
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FIG. 1: Shape of the functions: (a) secure return Ri a (xj), (b) 
non-secure return Ri n (xj), and (c) synergic function Si(xi). 



where ji > is a reactiveness to the appeal. When the 
synergic instinct is considered, the individual's reaction 
to the partner's appeal can be enhanced by love. For 
example, mothers often have a biased view of the beauty 
of their children. The synergic instinct [j| is described 
by 



V = {1 

Si(xi) = 



a t x. 




Si (xi)} jiAj, 
8 /(l + Xi 8 ) for x. 



(7) 



> 0, 
for Xi < 0, 



where the synergic function Si(xi) is an increasing and 
bounded function for Xi > 0, as shown in Fig. [TJc). 

Now, we can take the essential step for describing our 
love dynamical model. How does an individual know the 
partner's romantic feeling? In a real situation, the ro- 
mantic interaction is mediated by communication, e.g., a 
talk, a phone call, an email, a letter, etc. That is, time is 
required for the romantic feelings of someone to transfer 
to the other. In addition, Ackerman et al. very recently 
showed that an individual could delay one's confessing 
love to adjust potential costs and benefits [261 ] . 

Therefore, the oblivion, the return, and the instinct in 
the model of Eq. ([T]) do not proceed simultaneously, and 
the delay time r must be added to the return function. 
As a result, the love dynamical model can be described 
by DDEs as follows: 



ii(t) = Oi(xi(t)) + Ri(xj(t - r)) + Ufa®). (8) 

For the purpose of simplification, we consider the same 
delay time t for both individuals. 

It seems appropriate to comment on the recent work 
of Liao and Ran [25( . They introduced the time delay on 
return functions and observed the occurrence of a Hopf 
bifurcation. Using the same notation as Ref. (25|, their 
model is represented by 



±i(t) = -ciiXi{t) + bif(xj(t - Tj)) + jiAj 



(9) 



It contains the same basic processes of love dynamics: 
oblivion, return, and instinct. However, their return 
function bif(xj(t—Tj)) could not frilly consider the linear, 
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the secure, and the non-secure return functions. Also, 
they did not deal with the synergic instinct. As a result, 
they could not clarify in which type of return functions 
the Hopf bifurcation can arise. 

In their model, the condition for a Hopf bifurcation 
results in 



> — aia2 > b\b2 



„. dx 2 



(10) 



where (x*,x 2 ) is a fixed point of the model in Eq. (|9|). 
They also showed a numerical example of a Hopf bifur- 
cation in which ai = a.2 = 1, b\ = 1.5, b 2 = —2, and 
f(x) = tanh(z). Note that tanh(cc) is an increasing and 
bounded function, so individuals 1 and 2 exhibit a secure 
(&i > 0) and an anti-secure return (62 < 0), respectively. 
But an anti-secure individual seems unusual. 

Therefore, we will fully investigate the effect of time 
delay on the love dynamical models suggested by Rinaldi, 
Gragnani, and Feichtinger In the following section, 

we will explicitly show that, in both cases of synergic and 
non-synergic couples, if at least one individual exhibits 
a non-secure return, then the time delay on the return 
function can cause a Hopf bifurcation. 



When r = 0, all roots of Eq. (fT3|) have real negative 
parts if and only if the conditions 



(Hi) ai + a 2 > 0, ol\ol 2 — c\c 2 > 

hold. 

Now, let us assume that iuj (real positive w) is a root 
of Eq. (O. Then, we have 



CO — Ct\Ct 2 

(ai + a 2 )u} 
which lead to 



-C1C2 cos 2ojt, 
-C1C2 sin 2u>t, 



(14) 



p 2 + (ai 2 + a 2 2 )p + ai 2 a 2 2 - c x 2 c 2 2 = 0, (15) 
where p = ui 2 . It follows that if the conditions 



(H 2 ) a.\ + a 2 > 0, ai 2 a 2 - ci c 2 > 

are satisfied, then Eq. ([TBI has no positive roots. Hence, 
all roots of Eq. f| 13[) have real negative parts when r £ 
[0,oo). 

On the other hand, if the condition 



III. HOPF BIFURCATION ANALYSIS 

A. Non-synergic Couple 

First, let us consider a couple composed of non- 
synergic individuals. We follow the same steps in Refs. 
|13l [24j for verifying the occurrence of a Hopf bifurca- 
tion. The love dynamics for the non-synergic couple is 
described by 



x\(t) = -a 1 x 1 {t)+R 1 (x 2 {t-r))+'YiA 3 , (11) 
x 2 {t) = -a 2 x 2 (t) + R 2 {xi{t - t)) + 72 Ai. 

Here, we do not need to fix the type of return function 
Ri at this stage. It may be one of the linear, secure, and 
non-secure returns. We assume that {x\,x 2 ) is a fixed 
point of the model in Eq. (fTTj) . which is located in the 
first quadrant. Then, the linearization of the model in 



Eq. Jirj at (3 



lJ J, 2. 



is given by 



8x\{t) = —a\8x\{t) + c^Sxzit — r), 
Sx 2 (t) = — a 2 5x 2 (t) + c 2 Sxi(t - t), 



(12) 



where ci = dR\/dx 2 \ x * and C2 = dR 2 /dxi\ x *. The char- 
acteristic equation of Eq. (jT^j) is described by 



A 2 + (a\ + a 2 )X + axa 2 — c\C 2 e 



-2\t 



(13) 



(H 3 ) ai 2 a 2 2 - Cl 2 c 2 2 <0 

holds, then Eq. (|15[) has a unique positive root po = ujq 2 . 
Substituting uj into Eq. (fl"4"|) . we have the results 



CJ 



1 
1 

2uj 



\/(ai 2 -a 2 2 ) 2 +4ci 2 c 2 2 - («i 2 + a 2 2 ) 
2 



a\a 2 — uj - 



C\C 2 



nn 



1/2 



(16) 



for n = 0, 1, 2, • • • . Then, let us investigate the sign of 
Re [dX/dr]. The differentiation of Eq. (|13l) with respect 
to r and the substitutions r = tq and A = iu>o lead to 



Re 



dX 
dr 



T = T ,U=U) 



■iuj 4 + 2a; 2 (ai 2 + a 2 2 ) 
A 2 + B 2 



, (17) 



where 



A = (ai + 02) cos 2woT"o — 2wo sin2wo T o + 2ciC2To, 
B = (ai + CK2) sin 2ujqTq + 2luq cos 2u>oTq. 

It is easily obtained that 



Re 



dX 
fa 



> 0. 



(18) 



As a result, we have the following theorem from Corol- 
lary 2.4 in Ref. QJ]: 
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Theorem 1 

1. If (Hi) and (H 2 ) ZioZd, then the fixed point (x\,x 2 ) 
is asymptotically stable for all r > 0. 

i?. // (Hi) and (H3) ZioZd, then the fixed point (x\,x 2 ) 
is asymptotically stable for r < tq and unstable for 
t > tq. Furthermore, the love dynamical model in 
Eq. (1111) undergoes a Hopf bifurcation at (x\,x\) 
when t = tq. 

Note that the linear return Ri(xj) and the secure return 
Ri S (xj) are always increasing functions with a > 0, in 
contrast to the non-secure return. Because ot\ and a 2 are 
positive, the conditions (Hi) and (H 3 ) result in 

> — OL\OL 2 > C\C2- (19) 

Therefore, if no one exhibits a non-secure return, then 
the existence of time delay cannot disturb a steady state 
(x\,x 2 ). However, if at least one of them has a non- 
secure return and the inequality in Eq. (|T9"]) is satisfied, 
then the time delay on the return function leads to a 
Hopf bifurcation and a cyclic love dynamics. 



(H 4 ) a\ + a s 2 > 0, a{a 2 - c x c 2 > 

are satisfied. Then, assuming that iuj (real positive u) is 
a root of Eq. (|2"2l , we obtain 

p 2 + {(al) 2 + (^) 2 } P +(alf(a s 2 f - Cl 2 c 2 2 = 0. (23) 

where p — uj 2 . In the same manner as in Section [ill A[ it 
follows that if the conditions 

(H 5 ) K) 2 + (^) 2 >0, (a s 1 ) 2 (a s 2 ) 2 ~c 1 2 c 2 2 >0 

hold, then all roots of Eq. (|22[) have real negative parts 
when t£ [0,oo). If the condition 

(H 6 ) {a{) 2 {a s 2 ) 2 - Cl 2 c 2 2 <0 

is satisfied, then Eq. (|2"3"|) has a unique positive root 
Po = ujq 2 . Accordingly, we have the results 



B. Synergic Couple 

Second, let us investigate the effect of time delay on 
a couple composed of synergic individuals. In this case, 
the love dynamics is given by 

x\(t) = -axxx{t) + R x (x 2 (t-T)) (20) 

+ {l + S l {x 1 {i))} lx A 2 , 

x 2 {t) = -a 2 x 2 (t) + R 2 (xi(t - t)) 

+ {l + S 2 (x 2 (t))}i 2 A 1 . 

Also, the type of return function Ri is not fixed at this 
stage. It may be one of the linear, secure, and non-secure 
returns. Let [x\,x 2 ) be a fixed point of the model in 
Eq. (|2U)) . which is located in the first quadrant. The 
linearization of the model in Eq. (|2"0")) at (x\ ,x 2 ) is given 
by 



1 

71 



{(al) 2 - (a s 2 ) 2 } +4ci 2 c 2 2 (24) 
1/2 



{K) 2 + K) 2 } 



2uj 



a\a 2 — luq 



C\C 2 



n tv 
w ' 



for n = 0,1,2,- ••. 

Now, the remaining step is to fix the sign of Re [dX/dr] 
We can obtain 



Re 



dX 
d^ 



T=T ,U>=U> a 



4aV + 2aV{(o4) 2 + (al) 2 } 
C 2 + D 2 



. (25) 



where 



C = (aj + a 2 ) cos2woTo — 2wo sin2woTo + 2C1C2T0, 
D = (al +0^2) sin 2uj q t + 2wq cos 2uj t . 



Sxx(t) = -ax5xi(t) + c x 8x 2 (t - t) + diSxi(t%21) 
Sx 2 (t) = -a 2 5x 2 (t) + c 2 Sxi(t - t) + d 2 Sx 2 (t), 

where di = dSi/dxi\ x * ■ 71A2, d 2 = dS 2 /dx 2 \ x - ■ j 2 Ai. 
Then, the characteristic equation of Eq. (f2"Tj) is described 
by 



A 2 + (al + a 2 )X + a{a 2 - c\c 2 e 



-2\t 



0. 



(22) 



where af = on — di. For the case of r = 0, all roots of Eq. 
have real negative parts if and only if the conditions 



This leads to 



Re 



\dX 



> 0. 



T = T ,Ll) = U) 



(26) 



As a result, we have the following theorem from Corol- 
lary 2.4 in Ref. [lj]: 

Theorem 2 

1. If (H 4 ) and (H 5 ) hold, then the fixed point (x\,x 2 ) 
is asymptotically stable for all r > 0. 



5 



2. If (H4) and (He) hold, then the fixed point (x*, X2) 
is asymptotically stable for t < t and unstable for 
t > tq. Furthermore, the love dynamical model in 
Eq. (jlip undergoes a Hopf bifurcation at (x^x^) 
when t = tq. 

Note that the conditions (H4) and (He) result in 

> -Ka|| > cic 2 . (27) 

Thus, if no one exhibits a non-secure return, then the 
time delay cannot destabilize a steady state (x^x^)- On 
the other hand, if at least one of them has a non-secure 
return and the inequality in Eq. (|27j) is satisfied, then 
the time delay on the return function leads to a Hopf 
bifurcation and a cyclic love dynamics. Therefore, we 
obtain the same results as in the non-synergic case. 

IV. NUMERICAL BIFURCATION ANALYSIS 

In the previous section, we have theoretically proven 
the occurrence of a Hopf bifurcation. In this section, we 
show that a numerical bifurcation analysis supports the 
theoretical results on Hopf bifurcation and investigate 
additional bifurcation phenomena. The reactiveness to 
the love, i.e., ft, and the delay time r are considered as 
varying parameters, and the others are fixed. For numer- 
ical detection and continuation of a bifurcation poi nt in 
parameter space (ft, r), we use DDE-BIFTOOL [27| and 
KNUT [H- 

A. Non-synergic Couple 

First, let us investigate a non-synergic couple. Among 
various combinations, a couple composed of secure and 
non-secure individuals are considered as a proper case for 
observing the Hopf bifurcation based on the theoretical 
results. Thus, the love dynamics is described by 

x\(t) = -a 1 x 1 (t)+R 1 s (x 2 (t-T)) + 7 iA 2 , (28) 
a? a (t) = -a 2 x 2 (t)+R 2 n (xi(t-T)) +72^1, 

where the functional forms of the secure return R\ s and 
the non-secure return R 2 n are denoted in Eqs. (j4]) and 
([5]), respectively. 

Figure^a) is a bifurcation diagram in (ft, t). It shows 
a supercritical Hopf (sH), a limit point bifurcation of 
cycles (LPC) and two period doubling (PD) bifurcation 
curves. In this figure, the two PD curves represent the 
bifurcation from period-1 to period-2 for two different 
branches of the limit cycle. The black PD curve cor- 
responds to the period doubling of the stable limit cy- 
cle branch emerging from the LPC curve. On the other 
hand, the red one is associated with the period doubling 



of the limit cycle branch arising from the sH curve. Fig- 
ure EJa) also shows a chaotic region represented by gray 
dots. That is determined by the largest Lyapunov expo- 
nent Ai of the model in Eq. (|2"5)l . For this computation, 
we firstly obtain the time series of DDEs from the con- 
stant history function X\{t) = x 2 (t) = for — r < t < 0. 
Remember that we assume the two individuals to be com- 
pletely indifferent to each other when they first meet. 
Then, the largest Lyapunov exponent of the time series 
is estimated with the help of TISEAN [29J. The result 
of Fig. [2][a) also shows bistable phenomena: Two limit 
cycles emerging from the sH and the LPC curves coexist. 
Moreover, a limit cycle and a chaotic attractor coexist 
around the period-doubling bifurcation point '2' in Fig. 
EJa). It shows that the red PD curve is on the chaotic 
region. 

In Fig. [2jb), we plot the orbit diagram, i.e., the lo- 
cal maxima of the x\ variable as varying r at a fixed 
ft = 7.2. Here, the black and the red dots are obtained 
from xx{t) = x 2 (t) = and iKi(i) = 2A,x 2 (t) = 0.2 for 
— r < t < 0, respectively. The orbit diagram corresponds 
to the line indicated by blue arrows in Fig. [2ja) . In the 
following bifurcation diagrams, lines indicated by blue ar- 
rows correspond to their matching orbit diagrams. The 
result of Fig. [2]Jb) shows a cascade of period-doubling 
bifurcations. For a clearer illustration, we plot together 
the orbit diagram based on Poincare section (x 2 = 0) 
at the same fixed value ft = 7.2. It explicitly shows a 
period-doubling route to chaos. However, in the follow- 
ing orbit diagrams, we only plot those consisting of local 
extrema, because they are clearer than those based on 
Poincare section to understand how the bifurcation oc- 
curs. Also, Fig. EJb) and Fig. [2jc) agree well with the 
result that the time delay can induce a period doubling 
route to chaos 30]. 

It seems appropriate to comment on more PD curves. 
In Fig. [2ja), we could not present more PD curves corre- 
sponding to the bifurcation from period-2 to period-4 or 
for higher period because numerical detection and con- 
tinuation of a bifurcation point in DDEs is more subtle 
than that of ODEs. However, the result of the orbit 
diagram in Fig. Hfb) clearly explains that there exists 
a cascade of period doubling bifurcations in Fig. [2ja). 
Concerning this diagram, the non-zero history function 
does not agree with our real experience and the assump- 
tion on initial indifference. However, we add it to clarify 
the bistable phenomena observed in Fig. [2ja). The orbit 
diagram explicitly shows two coexisting limit cycles and 
the coexistence of a limit cycle and a chaotic attractor. 
The points '2' and '3' in Figs. Ufa) and[U[b) represent 
the same points. 

Now, let us show that the Hopf bifurcation in Fig. 
EJa) is supercritical. Though, for determining the di- 
rection of Hopf bifurcation, rigorous analysis based on 
the normal form method and the center manifold the- 
ory presented in Refs. [l3|, EJ is required, numeri- 
cal bifurcation analysis can be used to investigate it. 
In Fig. [3J we obtain the variation of the amplitude 
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FIG. 2: (a) Bifurcation diagram in (/3i,t) for the parameter 
value: ct\ = a.2 — 1, P2 = 1, and 71A2 = 72j4i = 0.5. (b) 
Orbit diagram of local maxima of x\ at a fixed value /3i = 7.2. 
(c) Orbit diagram based on Poincare section (X2 = 0). 
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FIG. 4: (a,b) Time series and (c) phase plot of the model in 
Eq. (|28[) f° r the parameter value: qi = 02 = 1, Pi = 7.2, 
/?2 = 1, 7i^4.2 = 72^4.1 = 0.5, and r = 3.2. The results are 
obtained from x\(t) = X2(t) = for t < 0. 




FIG. 3: Branch of the limit cycle emanating from the Hopf 
point and the steady state for the parameter value: a?i = 
«2 = 1, Pi = 7.2, p 2 = 1, and 71A2 = 72A1 = 0.5. The 
points '1' and '2' represent the same points in Figs. [2ja) and 

Hb). 



FIG. 5: (a) Bifurcation diagram in (/?2,t) for the parameter 
value: oei — a.2 — 1, Pi = 1, and 71A2 = 72^1 = 0.5. (b) 
Orbit diagram of local minima of X2 at a fixed value P2 = 7.2. 



Aiei = max(a;i(t))— min(j;i(t)) for the limit cycle branch 
emerging from the Hopf bifurcation point T' with in- 
creasing t at a fixed (3i = 7.2. The solid and the dashed 
curves represent the stable and the unstable branches, 
respectively. The result shows that a smooth transition 
from the steady state to the limit cycle arises across the 
Hopf point, supporting the Hopf bifurcation being su- 
percritical. The limit cycle branch arising from the Hopf 
point loses its stability at the period-doubling bifurcation 
point '2'. 

In Fig. |4j we show the time series and the phase plot of 
the model in Eq. (|2"5)) when it exhibits chaotic behavior 
(Ai ~ 0.038). The results tell us that the love dynamics 
of two individuals exhibits stormy patterns of feelings 



and a long-time unpredictable state. 

For the parameter space (02, t), we show the bifur- 
cation diagram in Fig. 02a). It consists of sH and PD 
curves including a chaotic region. In Fig. EJb), we plot 
the orbit diagram of local minima of Xi for fixed /?2 = 7.2. 
Here, the orbit diagram is obtained from 27 = X2 = for 
t < and clearly shows a period-doubling route to chaos. 
The results of Figs. EJa) and \5\b) are same as those of 
Figs. [Ha) and[UJb) except for the bistable phenomena. 
In both parameter spaces (/?i,t) and (^tt), we observe 
the same bifurcation structure consisting of a supercrit- 
ical Hopf and a cascade of period-doubling bifurcations, 
that is, a period-doubling route to chaos. 

Through a numerical bifurcation analysis, we have con- 
firmed the theoretical predictions on the occurrence of a 
Hopf bifurcation and obtained the following results: For 
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FIG. 6: (a) Bifurcation diagram in (/3i,t) for the parameter 
value: ai = ct2 = 1, /?2 = 1, 71^2 = 72^.1 = 0.5, and 
(Ji = <T2 = 0.5. (b) Orbit diagram of local minima of X2 at a 
fixed value /3i = 9.05. 
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FIG. 7: (a) Bifurcation diagram in (/?2,t) for the parameter 
value: Qi = «2 = 1, /?i = 1, 71A2 = 72 -A 1 = 0.5, and 
ai = (72 = 0.5. (b) Orbit diagram of local minima of X2 at a 
fixed value P2 = 4.3. 



ft smaller than ft c , a steady state of love dynamics can- 
not be destabilized by a Hopf bifurcation, no matter how 
long t is. That is, the existence of time delay does not 
influence a plateau of the love affair for insensitive cou- 
ple. ft c is a minimum value of ft, which satisfies the 
inequality of Eq. (|19[) for given parameter values. In 
those cases, ft c ~ 1.27 and ft> c ~ 1.68. For ft > ft c , 
the parameter value of r, where the supercritical Hopf 
bifurcation arises, becomes shorter as ft increases. For 
large ft and long t, the limit cycle undergoes a cascade of 
period doubling bifurcations resulting in chaotic motion, 
i.e., stormy patterns of feelings. 

B. Synergic Couple 

Let us now investigate the synergic couple. For the 
same reason as in the non-synergic case, the couple com- 
posed of a secure and a non-secure individuals are con- 
sidered. Then, the love dynamics is given by 

x\(t) = -aixi(t) + R 1 s (x 2 (t-T)) (29) 

+ {l + Si(xi{t))}~/iA 2 , 

x 2 (t) = -a 2 x 2 (t) + R 2 n (x 1 (t-T)) 

+ {l + S 2 (x 2 {t))}^ 2 A u 

where the functional forms of the secure return R\ s and 
the non-secure return R 2 n are described by Eqs. ((U) and 
([5]), respectively. Here, the synergic functions Si and S 2 
are given by Eq. (J7J . 

For the parameter space (ft , r) , we plot the bifurca- 
tion and the orbit diagrams of the model in Eq. (f2U| in 
Figs. [H[a) and[5Jb), respectively. In Fig. [SJb), the black 
and the red dots are obtained from xi(t) = x 2 (t) = 
and xi(t) = 2A,x 2 (t) = 0.2 for t < 0, respectively. The 
above results exhibit bistable phenomena. The point '1' 
in Figs. [Hla) andJBJb) represents the same point, which 
corresponds to a period doubling of the limit cycle branch 



arising from the sH curve. However, we cannot continue 
other bifurcation curves involved in the bistable phenom- 
ena, such as the LPC and the black PD curves in Fig. 
Ha). 

In Figs. [Ha) andJTJb), we present the results of the 
numerical bifurcation analysis for the parameter space 
W2,t). Here, the orbit diagram is obtained from xi(t) = 
x 2 (t) = for t < 0. Equivalently to the non-synergic 
case, the results of Figs. EJa) and [TJb) are the same as 
those for Figs. [6ja) andHJb) except for the bistable phe- 
nomena. Therefore, we explicitly show that a supercrit- 
ical Hopf and a cascade of period-doubling bifurcations, 
i.e., a period-doubling route to chaos, is a universal bi- 
furcation structure in our models in Eqs. (|28|) and (|29|) . 
For the synergic case, we identify that ft c ~ 0.79 and 
ft c ~ 0.60. 

V. CONCLUSION 

In conclusion, we have investigated the effect of time 
delay on a simplified mathematical model that describes 
the dynamics of love between two individuals in the field 
of social science. By analyzing the characteristic equa- 
tion of linearization of the model, we have theoretically 
proven that if no one exhibits a non-secure return in both 
cases of synergic and non-synergic couples, then the exis- 
tence of time delay cannot disturb a steady state of love 
dynamics. However, if at least one of them has a non- 
secure return, then the time delay on the return function 
can cause a Hopf bifurcation and a cyclic love dynam- 
ics. Through a numerical bifurcation analysis, we have 
confirmed the theoretical predictions on the occurrence 
of the Hopf bifurcation and obtained the universal bifur- 
cation structure consisting of a supercritical Hopf and 
a cascade of period-doubling bifurcations, resulting in 
chaotic motion, which exhibits stormy patterns of feel- 
ings and a long-time unpredictable state. We have also 
ascertained that the existence of time delay does not in- 
fluence a steady state of the love affair for an insensitive 
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couple. 

Through further investigation, we hope that the time 
delay effect on a more realistic love dynamical model, 
including a love triangle, can be investigated. Moreover, 
it seems that such modeling approach with time delay 
can be applied to various dynamical phenomena in the 
field of social science, including a study of the proper 
relation between supplier and consumer in the field of 
management. 
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